{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HW3: Problem 3\n",
    "\n",
    "The code here is used for problem 3 of HW3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import statsmodels.formula.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data\n",
    "df = pd.read_csv('../data/Data_Lecture_3.csv', index_col=0, parse_dates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>AAL</th>\n",
       "      <th>AAP</th>\n",
       "      <th>AAPL</th>\n",
       "      <th>ABC</th>\n",
       "      <th>ABT</th>\n",
       "      <th>ACN</th>\n",
       "      <th>ADBE</th>\n",
       "      <th>ADI</th>\n",
       "      <th>ADM</th>\n",
       "      <th>...</th>\n",
       "      <th>XYL</th>\n",
       "      <th>YUM</th>\n",
       "      <th>ZBH</th>\n",
       "      <th>ZION</th>\n",
       "      <th>MKT</th>\n",
       "      <th>SMB</th>\n",
       "      <th>HML</th>\n",
       "      <th>RF</th>\n",
       "      <th>CPI</th>\n",
       "      <th>IndProd</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2013-12-31</th>\n",
       "      <td>216.0</td>\n",
       "      <td>423.0</td>\n",
       "      <td>293.0</td>\n",
       "      <td>63.0</td>\n",
       "      <td>407.0</td>\n",
       "      <td>150.0</td>\n",
       "      <td>124.0</td>\n",
       "      <td>349.0</td>\n",
       "      <td>116.0</td>\n",
       "      <td>343.0</td>\n",
       "      <td>...</td>\n",
       "      <td>201.0</td>\n",
       "      <td>126.0</td>\n",
       "      <td>277.0</td>\n",
       "      <td>267.0</td>\n",
       "      <td>218.0</td>\n",
       "      <td>59.0</td>\n",
       "      <td>42.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>36.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-01-31</th>\n",
       "      <td>238.0</td>\n",
       "      <td>430.0</td>\n",
       "      <td>375.0</td>\n",
       "      <td>214.0</td>\n",
       "      <td>401.0</td>\n",
       "      <td>119.0</td>\n",
       "      <td>138.0</td>\n",
       "      <td>405.0</td>\n",
       "      <td>165.0</td>\n",
       "      <td>390.0</td>\n",
       "      <td>...</td>\n",
       "      <td>210.0</td>\n",
       "      <td>147.0</td>\n",
       "      <td>212.0</td>\n",
       "      <td>240.0</td>\n",
       "      <td>234.0</td>\n",
       "      <td>73.0</td>\n",
       "      <td>47.0</td>\n",
       "      <td>46.0</td>\n",
       "      <td>53.0</td>\n",
       "      <td>54.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-02-28</th>\n",
       "      <td>368.0</td>\n",
       "      <td>452.0</td>\n",
       "      <td>399.0</td>\n",
       "      <td>179.0</td>\n",
       "      <td>379.0</td>\n",
       "      <td>129.0</td>\n",
       "      <td>128.0</td>\n",
       "      <td>396.0</td>\n",
       "      <td>121.0</td>\n",
       "      <td>257.0</td>\n",
       "      <td>...</td>\n",
       "      <td>236.0</td>\n",
       "      <td>84.0</td>\n",
       "      <td>267.0</td>\n",
       "      <td>210.0</td>\n",
       "      <td>237.0</td>\n",
       "      <td>98.0</td>\n",
       "      <td>41.0</td>\n",
       "      <td>48.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>68.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-03-31</th>\n",
       "      <td>337.0</td>\n",
       "      <td>448.0</td>\n",
       "      <td>396.0</td>\n",
       "      <td>218.0</td>\n",
       "      <td>309.0</td>\n",
       "      <td>162.0</td>\n",
       "      <td>142.0</td>\n",
       "      <td>403.0</td>\n",
       "      <td>144.0</td>\n",
       "      <td>227.0</td>\n",
       "      <td>...</td>\n",
       "      <td>372.0</td>\n",
       "      <td>85.0</td>\n",
       "      <td>251.0</td>\n",
       "      <td>250.0</td>\n",
       "      <td>237.0</td>\n",
       "      <td>97.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>49.0</td>\n",
       "      <td>61.0</td>\n",
       "      <td>68.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-04-30</th>\n",
       "      <td>334.0</td>\n",
       "      <td>450.0</td>\n",
       "      <td>393.0</td>\n",
       "      <td>238.0</td>\n",
       "      <td>226.0</td>\n",
       "      <td>104.0</td>\n",
       "      <td>55.0</td>\n",
       "      <td>379.0</td>\n",
       "      <td>239.0</td>\n",
       "      <td>289.0</td>\n",
       "      <td>...</td>\n",
       "      <td>314.0</td>\n",
       "      <td>151.0</td>\n",
       "      <td>244.0</td>\n",
       "      <td>257.0</td>\n",
       "      <td>223.0</td>\n",
       "      <td>110.0</td>\n",
       "      <td>81.0</td>\n",
       "      <td>62.0</td>\n",
       "      <td>73.0</td>\n",
       "      <td>82.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 455 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                A    AAL    AAP   AAPL    ABC    ABT    ACN   ADBE    ADI  \\\n",
       "Date                                                                        \n",
       "2013-12-31  216.0  423.0  293.0   63.0  407.0  150.0  124.0  349.0  116.0   \n",
       "2014-01-31  238.0  430.0  375.0  214.0  401.0  119.0  138.0  405.0  165.0   \n",
       "2014-02-28  368.0  452.0  399.0  179.0  379.0  129.0  128.0  396.0  121.0   \n",
       "2014-03-31  337.0  448.0  396.0  218.0  309.0  162.0  142.0  403.0  144.0   \n",
       "2014-04-30  334.0  450.0  393.0  238.0  226.0  104.0   55.0  379.0  239.0   \n",
       "\n",
       "              ADM  ...    XYL    YUM    ZBH   ZION    MKT    SMB   HML    RF  \\\n",
       "Date               ...                                                         \n",
       "2013-12-31  343.0  ...  201.0  126.0  277.0  267.0  218.0   59.0  42.0  32.0   \n",
       "2014-01-31  390.0  ...  210.0  147.0  212.0  240.0  234.0   73.0  47.0  46.0   \n",
       "2014-02-28  257.0  ...  236.0   84.0  267.0  210.0  237.0   98.0  41.0  48.0   \n",
       "2014-03-31  227.0  ...  372.0   85.0  251.0  250.0  237.0   97.0  43.0  49.0   \n",
       "2014-04-30  289.0  ...  314.0  151.0  244.0  257.0  223.0  110.0  81.0  62.0   \n",
       "\n",
       "             CPI  IndProd  \n",
       "Date                       \n",
       "2013-12-31  34.0     36.0  \n",
       "2014-01-31  53.0     54.0  \n",
       "2014-02-28  60.0     68.0  \n",
       "2014-03-31  61.0     68.0  \n",
       "2014-04-30  73.0     82.0  \n",
       "\n",
       "[5 rows x 455 columns]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = df.rolling(11).sum().shift()\n",
    "# Drop NaN values (first 11 months)\n",
    "df = df.dropna()\n",
    "# Rank in ascending order, for each time step\n",
    "ranked_df = df.rank(axis=1)\n",
    "ranked_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Date\n",
       " 2013-12-31    1.239963\n",
       " 2014-01-31   -1.529257\n",
       " 2014-02-28   -2.585661\n",
       " 2014-03-31   -2.791504\n",
       " 2014-04-30   -4.220694\n",
       " dtype: float64, Date\n",
       " 2013-12-31    55.403097\n",
       " 2014-01-31    50.718519\n",
       " 2014-02-28    49.152139\n",
       " 2014-03-31    51.018520\n",
       " 2014-04-30    47.913327\n",
       " dtype: float64)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Find the 100 companies with best and worst returns\n",
    "worst_mask = ranked_df<=100\n",
    "best_mask = (355<ranked_df) & (ranked_df<=455)\n",
    "worst_avg = df.where(worst_mask).mean(axis=1)\n",
    "best_avg = df.where(best_mask).mean(axis=1)\n",
    "worst_avg.head(), best_avg.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Date\n",
       "2013-12-31    54.163134\n",
       "2014-01-31    52.247775\n",
       "2014-02-28    51.737800\n",
       "2014-03-31    53.810023\n",
       "2014-04-30    52.134021\n",
       "dtype: float64"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get the difference between best and worst averages\n",
    "factor = best_avg - worst_avg\n",
    "factor.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAfK0lEQVR4nO3dfbRcdX3v8feHJOCBigckYDiYgi033FYuSTiyanMXVVARdAFi64XbIiLe2FWx2EoweG+tyvISwVatXZc2Pl1WfaiKIbpEeTABb4XK8oQEQSFFHlQOkQQkhkoqIXzvH3ufyeQwZ57Ofpz5vNaadebsmdnz/c2cs7/797gVEZiZmQHsU3YAZmZWHU4KZmbW4KRgZmYNTgpmZtbgpGBmZg1zyw6gG4ccckgceeSRZYdhZlYrGzZseCwi5vfymlokhSOPPJKJiYmywzAzqxVJP+n1NW4+MjOzBicFMzNrcFIwM7MGJwUzM2twUjAzs4ZajD4yMxtkazdOcuUNm3lk+04OHx1hxSmLOHPJWCmxOCmYmZVo7cZJLl1zFzt37QZgcvtOLl1zF0ApicHNR2ZmJbryhs2NhDBl567dXHnD5lLicVIwMyvRI9t39rQ9b7klBUmLJG1quu2Q9C5J75c02bT9tLxiMDOrusNHR3ranrfckkJEbI6IxRGxGDgeeAq4Nn34o1OPRcQ384rBzKzqVpyyiJF5c/baNjJvDitOWVRKPEV1NJ8M3B8RP5FU0FuamVXfVGfysI0+Ohv4YtPvF0p6MzABvDsinpj+AknLgeUACxcuLCRIM7MynLlkrLQkMF3uHc2S9gVOB76SbroK+C1gMbAF+JtWr4uI1RExHhHj8+f3tPKrmZn1qYjRR6cCd0TEowAR8WhE7I6IZ4FPAicUEIOZmXWhiKRwDk1NR5IWND32BuDuAmIwM7Mu5NqnIGl/4NXA25s2XyFpMRDAQ9MeMzOzEuWaFCLiKeCF07adm+d7mplZ/zyj2czMGpwUzMyswUnBzMwanBTMzKzBScHMzBqcFMzMrMFJwczMGpwUzMyswUnBzMwanBTMzKzBScHMzBqcFMzMrMFJwczMGpwUzMyswUnBzMwanBTMzKzBScHMzBqcFMzMrCG3pCBpkaRNTbcdkt4l6WBJN0m6L/15UF4xmJlZb3JLChGxOSIWR8Ri4HjgKeBaYCWwLiKOBtalv5uZWQUU1Xx0MnB/RPwEOAO4Ot1+NXBmQTGYmVkHRSWFs4EvpvcPi4gtAOnPQ1u9QNJySROSJrZt21ZQmGZmwy33pCBpX+B04Cu9vC4iVkfEeESMz58/P5/gzMxsL0XUFE4F7oiIR9PfH5W0ACD9ubWAGMzMrAtFJIVz2NN0BPB14Lz0/nnA1wqIwczMupBrUpC0P/BqYE3T5lXAqyXdlz62Ks8YzMyse3Pz3HlEPAW8cNq2x0lGI5mZWcV4RrOZmTU4KZiZWYOTgpmZNTgpmJlZg5OCmZk1OCmYmVmDk4KZmTU4KZiZWYOTgpmZNTgpmJlZg5OCmZk1OCmYmVmDk4KZmTU4KZiZWUOuS2ebmQ2btRsnufKGzTyyfSeHj46w4pRFnLlkrOywuuakYGaWkbUbJ7l0zV3s3LUbgMntO7l0zV0AtUkMbj4yM8vIlTdsbiSEKTt37ebKGzaXFFHv8r4c56ikayTdK+keSS+X9H5Jk5I2pbfT8ozBzKwoj2zf2dP2Ksq7pvBx4PqIOAY4Drgn3f7RiFic3r6ZcwxmZoU4fHSkp+1VlFtSkHQgcCLwaYCIeDoituf1fmZmZVtxyiJG5s3Za9vIvDmsOGVRSRH1Ls+awkuAbcBnJW2U9ClJB6SPXSjpB5I+I+mgVi+WtFzShKSJbdu25RimmVk2zlwyxuVnHcvY6AgCxkZHuPysY2vTyQygiMhnx9I48D1gWUTcLunjwA7g74HHgAAuAxZExFvb7Wt8fDwmJiZyidPMbFBJ2hAR4728Js+awsPAwxFxe/r7NcDSiHg0InZHxLPAJ4ETcozBzMx6kFtSiIifAz+TNNWYdjLwI0kLmp72BuDuvGIwM7Pe5D157Z3A5yXtCzwAnA/8naTFJM1HDwFvzzkGMzPrUq5JISI2AdPbs87N8z3NzKx/ntFsZmYNXvvIzKyCylpYz0nBzKxiylxYz81HZmYVU+bCek4KZmYVU+bCek4KZmYVU+bCek4KZmYVU+bCeu5oNjOrmKnOZI8+MjMzIEkMZayu2jYpSNoH+MOI+HJB8ZhZSep+wfl+DGOZO2nbp5CuZHphQbGYWUmmxsVPbt9JsGdc/NqNk2WHlpthLHM3uulovknSxZJeLOngqVvukZlZYQbhgvO9GsYyd6ObPoWpC+C8o2lbkFxZzcwGwCBccL5Xw1jmbnRMChFxVBGBWH24HXbwHD46wmSLg2GdLjjfq2Esczc6Nh9JmifpzyVdk94ulDSviOCseopsh127cZJlq9Zz1MrrWLZq/dC39U7J43MZhAvO92oYy9yNbpqPrgLmAf8n/f3cdNvb8grKqqtdO2yWtYU8FwSrc00nr8+lzHHxZRnGMnejm6Twsog4run39ZLuzCsgq7ai2mHzSj5lrj6ZhTyTclnj4ss0jGXupJuksFvSb0XE/QCSXgLs7vAaG1BFtcPmlXyKqunkpcjO0aJqVHWuuQ2iboakrgBulnSLpO8A64GLu9m5pNG0H+JeSfdIenk6pPUmSfelPw+aTQGsWEW1w+a1IFjdR5wUtVBaUX1HnitQPd0khe8CRwN/nt4WAbd2uf+PA9dHxDHAccA9wEpgXUQcDaxLf7eaOHPJGJefdSxjoyMIGBsd4fKzjs38zC6v5FPm6pNZKCopFzWGfxjmCtRtwEQ3zUf/GhFLgR9MbZB0B7C03YskHQicCLwFICKeBp6WdAbwivRpVwO3AO/pMW4rURHtsHl1Aq44ZdFefQpQrxEnRXWOFlWjqnvNrZM69mHNmBQkvQgYA0YkLQGUPnQgsH8X+34JsA34rKTjgA3ARcBhEbEFICK2SDp0hvdfDiwHWLhwYXelsYGSR/IZhBEnRSTlovqOBn2uQB37sNrVFE4hOcs/Avjbpu07gPd2ue+lwDsj4nZJH6eHpqKIWA2sBhgfH49uX2fVUOXOQ4846ayoGlXda26d1LEmNGNSiIirgaslvTEivtrHvh8GHo6I29PfryFJCo9KWpDWEhYAW/vYt1VYHavMtreialSDUHNrp441IUW0PwmX9L+BKyJie/r7QcC7I+J/ddy59C/A2yJis6T3AwekDz0eEaskrQQOjohL2u1nfHw8JiYmOpfGKmHZqvUt/xHGRke4deVJJURkVa65DbLpJ0iQ1ITeePwYN9+7LffvQ9KGiBjv5TXddDSfGhGN5qKIeELSaUDHpAC8E/i8pH2BB4DzSUY8fVnSBcBPgT/qJWCrvkEcS19nrrmVp1VN6JXHzOerGyYr+310kxTmSNovIn4NIGkE2K+bnUfEJqBVljq5+xCt6qYfmEf3n8cTT+16zvPyGktf1X+uqqhjZ+cgmd6HtWzV+kp/H90khc8B6yR9lmTJ7LeSDCU1a3lgnrePmDdH7Nq9p2kyi87D6cnnqaefqfQ/V1XUsbNzkFX9++hm6ewrJN1FcnYv4LKIuCH3yKwWWp2F7no2GB2ZxwH7zc2sWadV8plJVf65qqKOnZ2DrOrfRzc1BSLiW8C3co7FamimA/Avd+5i01+/JrP3aZV8ZtL8z+U+h8Ef9lk3Vf8+urmewu9J+r6kf5f0tKTdknYUEZxVX1HLRnR79t/8z+V1dRJFLU1i3an699FNTeHvgbOBr5B0Gr8Z+O08g7L6KOqsZ6Yqd7tmKnew7uEJe9VS5e+j2+ajH0uaExG7SZatuC3nuKwmipp8NFPyef/pvzvje1W9Q8+sirpJCk+l8ww2SboC2MKeSWhmmZ31tGv/7yf5VL1DbxCV1YfjvqPsdJMUziXpe7gQ+AvgxcAb8wzKhk83cw56TT5V79DLSlUOiGXNG/F8lWy1WyV1XUScDPxZRLwH+A/gA4VFZkMlj/b/QV9XB8o9IFZl3oj7jrLVrqawQNIfAKdL+mf2LJ0NQETckWtkNlTyav+vcodeFso6IFZp3oj7jrLVLim8j2RV0yOAv2HvpBCAVzazzLj9vz9lHRB7mTeyj8RRK6/Lrabmv51stVs6+xrgGkl/FRGXFRiTDYHpTQ/TFwmDwWz/z1pZB8Reks7udCXmvJq2iuw7qkr/TZ46Tl5zQrCstZpU9tUNk7zx+LHKTuipqqKu2TzdTElndGRe4zucIz3n8Tyuv1zUZLBhmQzZ1TwFsyzN1A5+873bfL2FHpXVmd7NvJGjVl7X8rVZNG21OmPv9Lcz27P8YenQdlKwwrljMFtldKZ3k4xmatp6wcg8lq1a3/fBuZ8RV1mM0hqWv9t2Q1IPbvfCiPhF9uHYMBjWjsEqtUdnEUunZNSqNjFvH/Grp59h+87kehv9HJz7OWPP4ix/WP5u2/UpPAZsAibS24amm6+NaX0rqx28TFVqjy4qllZt/b/xvLl7XWcDeu9n6OeMPYuz/GH5u23XfPQJ4BXArcAXge9Gpws6TyPpIeBJYDfwTESMp9dq/h/AtvRp742Ib/YWttXN9DPToq5RWxVVao8uMpbptYks+hn6OWPP4ix/GCZDQvshqRdJEkliOBf4hKQbgasi4sEe3uOVEfHYtG0fjYiP9Byt1VKr9tyvbpgcqtFFVWqPLjOWbg7OnZq2uh2C2ryfF4zMy+RqgIM+GRI6DEmNxM3AJcA/AOcDryoiMBsc7c5Mh0VR153oRpmxdGqC6aZpq5shqNP3s33nLgg4aP95HvLcQbuO5gOAM4D/BswH1gBLI+JnPew/gBslBfCPEbE63X6hpDeT9E28OyKe6Ct6q4UqnSWXpUqL85UZS6cmmG6btjqdsc90mdj9953LxvfNfEXAKg0GKEu7PoWtwH0k/Qk/JjnAv0zSywAiYk0X+18WEY9IOhS4SdK9wFXAZen+LiNZQuOt018oaTmwHGDhwoVdF8iqp+qjNoo4EBTdHp31MuRZandAz+oEop/9eLXVRLuk8BWSA/cx6a1ZkNQc2oqIR9KfWyVdC5wQEf9v6nFJnwS+McNrVwOrAcbHx3vq4AZn/Cqp0lnydEUeCIpqj85jGfKiZHUC0c9+qjQYoEwz9ilExFsi4vwZbs85s59O0gGSnj91H3gNcLekBU1PewNw92wLMV2Vhv9Zta9JO4j9HXUuU1bDPvvZz0y1iMntOzlq5XUsW7V+KI4hbWc0S3opsAL4XZLawY+Aj0TEXV3s+zDg2mQAE3OBL0TE9ZL+SdLidH8PAW/vP/zW+s34rl20VsREp7Lk2d9R1t9TnftwsmrayvJKfcBeJ5fN+x9E7TqazwA+AlzOnqWzjwfWSLo4Ir7WbscR8QBwXIvt584q4i64PbG1fg5Sdftcei1jXv0dZX5uVe/D6aTVCUQ/f7tZXKlvumFoTmo3JPWDwKsj4jMR8YOIuDMiPgO8On2ssvoZclfnKnc3+m1Sq9Pn0k8Z85qlWubnNmgzb8uagT2TOtS4ZqNdUpgXEQ9N35hum5dXQFnIsj1xUP4A+j1I1elz6aeMefV3lPm5VbkPpx9FJtgzl4xx68qTeHDV6xir0NySIrXrU9glaWFE/LR5o6TfBJ7JN6zZybI9cVD+APo9SNXpc+m3jHn0d5T9uVW1D6cfZSXYfmZOD0JfZLuawl8D35b0FknHSnqppPOBG0ku1VlpzRn/1pUndfySBq3KPV2/s1jr9LlUadZwnT63qivre+1n5vQgjHRst/bRWkkPAu8G3knS0fxD4E0RcWdB8RWm7Ak9eet3rkCdPpcqzYeo0+fWrbLOiMuegd3rzOm6d0arx4VPSzE+Ph4TE16te7YGrZrbyjCUsQzTR1NBcmAuqq+iqt/rUSuvo9URVMCDq15XdDjPjUPaEBHjPb1mpqQg6evtXhgRp/fyRrORVVKo6h9WlfgzKlZdPu9lq9a37CMZGx0Z6kuoVv1z6ScptOtofjnwM5K1j26HtqO0Kq/fS/jV4R82K3Wbk1B3dfq86zQKrUhVarLMSruO5hcB7wVeCnycZH7CYxHxnYj4ThHBZanXYW2D2IHUSZ3mJAyCOn3eVerEr5JBG/4L7TuadwPXA9dL2g84B7hF0gcj4hNFBZiVTuua9LuE7yDx2WCx6vR5D+IZcVYGafgvdLjIjqT9JJ0FfA54B/B3dLE6ahW1O6NpVROo0z9sVnw2WKw6fd6DeEZsrc2YFCRdDdwGLAU+EBEvi4jLIqKW7Setxo1P11x1r9M/bFY8tr5Ydfu8e537Y/XUrqZwLvCfgIuA2yTtSG9PStpRTHjZ6XVdk7r9w2bBZ4PF8udtVTS08xS6GUo2bKOPzGywZD0kdaB103E2aB1IZmadDG1SGMRlCMzMZmtokwK4JmBmNt1QJ4VB4z6QcmXx+fs7tLLlmhQkPQQ8CewGnomIcUkHA18CjiS5RvObIuKJPOMYBnVaMmEQZfH5+zu0Kmg7eS0jr4yIxU094CuBdRFxNLAu/d1mqU5LJgyiLD5/f4dWBUUkhenOAK5O718NnFlCDANnGGdgV0kWn7+/Q6uCvJNCADdK2iBpebrtsIjYApD+PLTVCyUtlzQhaWLbtm05h1l/wzgDu0qy+Pz9HVoV5J0UlkXEUuBU4B2STuz2hRGxOiLGI2J8/vz5+UU4IIZxBnaVZPH5+zu0Ksi1ozkiHkl/bpV0LXAC8KikBRGxRdICYGueMQwLz7soVxafv79Dq4LclrmQdACwT0Q8md6/CfggcDLweESskrQSODgiLmm3L1+O08ysd1Vb5uIw4FpJU+/zhYi4XtL3gS9LugD4KfBHOcZgZmY9yC0pRMQDwHEttj9OUlswM7OKKWNIqpmZVZSXuRhgXjLBzHrlpDCgvGTC7Dih2rByUhhQ7ZZM8MGtvSolVCcnK5r7FAaUl0zoX1XWIJpKTpPbdxLsSU5rN9byMulWE04KA8pLJvSvKgm1KsnJhouTwoDykgn9q0pCrUpysuHipDCgzlwyxuVnHcvY6AgCxkZHuPysY90e3YWqJNSqJCcbLu5oHmC+3Gh/qrIG0YpTFu3V4Q2u7Vn+nBTMWqhCQq1KcrLh4qRgVmFVSE42XJwUzLrg+QI2LJwUzDqo0mQ2s7x59JFZB54vYMPENYWacPNFeTxfwIaJk0IN5Nl84WTT2eGjI0y2SACeL2CDyM1Hs7R24yTLVq3nqJXXsWzV+lzWpcmr+cJr63SnKpPZzIrgpDALRR1U82q+cFt5d1rNDn/j8WNcecPmXE8GzMqQe/ORpDnABDAZEa+X9H+BPwB+mT7lLRGxKe848lDU8tR5NV+4rbx7zfMFPBrJBlkRNYWLgHumbVsREYvTWy0TAhR3UM2r+cJr6/THNSwbZLkmBUlHAK8DPpXn+5SlqINqXovbua28P65h2SDLu/noY8AlwPOnbf+QpPcB64CVEfHr6S+UtBxYDrBw4cKcw+xPkQuW5bHcgdfW6Y9HI9kgU0Tks2Pp9cBpEfFnkl4BXJz2KSwAfg7sC6wG7o+ID7bb1/j4eExMTOQS52x5SOfwmd6nAMnJgJcmt6qRtCEixnt5TZ41hWXA6ZJOA54HHCjpcxHxJ+njv5b0WeDiHGPInRcsGz6uYdkgy62msNebTKspRMQWSQI+CvxHRKxs9/oq1xTMzKqqajWFmXxe0nxAwCbgT0uIwczMWigkKUTELcAt6f2TinhPMzPrnWc0m5lZg5OCmZk1eJXUHnkIqpkNMieFHnjNGzMbdE4KPehmAbwiaxKutZhZ1pwUetBpzZsiaxKutZhZHtzR3INOC+AVuXqmV+o0szw4KfSg06qiRa6e6ZU6zSwPTgo96LSEdZHXJ/C1EMwsD+5T6FG7BfCKXEq7yPcys+HhpJChIlfP9EqdZpaHQlZJnS2vkmpm1ru6rJJqXfAcBDMrg5NCAXo9wHsOgpmVxaOPcjZ1gJ/cvpNgzwF+7cbJGV/jOQhmVhYnhZz1c4D3HAQzK4uTQs76OcB7DoKZlSX3pCBpjqSNkr6R/n6UpNsl3SfpS5L2zTuGMvVzgO80c9rMLC9F1BQuAu5p+v3DwEcj4mjgCeCCAmIoTT8H+E4zp83M8pLr6CNJRwCvAz4E/KUkAScB/z19ytXA+4Gr8oyjTP1OMms3c9rMLC95D0n9GHAJ8Pz09xcC2yPimfT3h4GWRz5Jy4HlAAsXLsw5zHz5AG9mdZFb85Gk1wNbI2JD8+YWT205pToiVkfEeESMz58/P5cYzcxsb3nWFJYBp0s6DXgecCBJzWFU0ty0tnAE8EiOMVSSZyubWVXlVlOIiEsj4oiIOBI4G1gfEX8M3Az8Yfq084Cv5RVDFfUzmc3MrChlzFN4D0mn849J+hg+XUIMpfFsZTOrskLWPoqIW4Bb0vsPACcU8b5VNNOktcntO1m2ar2blMysVJ7RXLCZJq0J3KRkZqVzUihYq8ls4rlDsNykZGZlcFIoWKvZyjNd5sgL4JlZ0Xw9hRJMn8y2bNV6JlskAC+AZ2ZFc02hArwAnplVhWsKFdDv+khmZllzUqgIr49kZlXg5iMzM2twUjAzswYnBTMza3BSMDOzBicFMzNrUMRM82mrQ9I24Cd9vvwQ4LEMwynboJUHBq9MLk+1DVN5fjMierpKWS2SwmxImoiI8bLjyMqglQcGr0wuT7W5PO25+cjMzBqcFMzMrGEYksLqsgPI2KCVBwavTC5Ptbk8bQx8n4KZmXVvGGoKZmbWJScFMzNrqF1SkPRiSTdLukfSDyVdlG4/WNJNku5Lfx6Ubj9G0r9K+rWki1vsb46kjZK+UXRZmmLIrEySHpJ0l6RNkiYGoDyjkq6RdG+6v5fXtTySFqXfy9Rth6R31bU86WN/ke7jbklflPS8mpfnorQsPyzju2mKo9cy/bGkH6S32yQd17Sv10raLOnHklZ2fPOIqNUNWAAsTe8/H/g34HeAK4CV6faVwIfT+4cCLwM+BFzcYn9/CXwB+MYglAl4CDhkUL4j4Grgben9fYHROpenaZ9zgJ+TTC6qZXmAMeBBYCT9/cvAW2pcnpcCdwP7k1xW4NvA0UWXp88y/T5wUHr/VOD2pr+z+4GXpP8/dwK/0+69a1dTiIgtEXFHev9J4B6SP84zSA4gpD/PTJ+zNSK+D+yavi9JRwCvAz5VQOgzyrJMVZBVeSQdCJwIfDp93tMRsb2QQjTJ6fs5Gbg/Ivqdqd+3jMszFxiRNJfkYPpIzuE/R4bl+c/A9yLiqYh4BvgO8IYCivAcfZTptoh4It3+PeCI9P4JwI8j4oGIeBr453QfM6pdUmgm6UhgCXA7cFhEbIHkAyU5G+jkY8AlwLM5hdizDMoUwI2SNkhanlec3ZpleV4CbAM+mzbxfUrSATmG21EG38+Us4EvZh1fr2ZTnoiYBD4C/BTYAvwyIm7MM95OZvn93A2cKOmFkvYHTgNenF+03emjTBcA30rvjwE/a3rs4XTbjGqbFCT9BvBV4F0RsaOP178e2BoRGzIPrk+zLVNqWUQsJalCvkPSiZkF2KMMyjMXWApcFRFLgF+RVJlLkdH3g6R9gdOBr2QVW59xzPZ/6CCSs86jgMOBAyT9SbZR9hTPrMoTEfcAHwZuAq4naWp5JtMge9RrmSS9kiQpvGdqU4untZ2HUMukIGkeyQf1+YhYk25+VNKC9PEFwNYOu1kGnC7pIZIq1UmSPpdTyB1lVCYi4pH051bgWpLqY+EyKs/DwMMRcXv6+zUkSaJwWX0/qVOBOyLi0ewj7U5G5XkV8GBEbIuIXcAakrbtwmX4//PpiFgaEScCvwDuyyvmTnotk6T/QtIUfkZEPJ5ufpi9aztH0KGJr3ZJQZJI2pjviYi/bXro68B56f3zgK+1209EXBoRR0TEkSRV+fURUcpZTlZlknSApOdP3QdeQ1IlLlSG39HPgZ9JWpRuOhn4UcbhdpRVeZqcQ4lNRxmW56fA70naP93nySRt34XK8vuRdGj6cyFwFiV9T72WKY13DXBuRPxb0/O/Dxwt6ai0hnp2uo+ZzbaXvOgb8F9Jqj8/ADalt9OAFwLrSDL7OuDg9PkvIsmWO4Dt6f0Dp+3zFZQ7+iiTMpG0wd+Z3n4I/M86lyd9bDEwke5rLekIixqXZ3/gceAFdf97Sx/7AHAvycnHPwH71bw8/0Jy4nEncHKNvqNPAU80PXeiaV+nkYxeur+bY4KXuTAzs4baNR+ZmVl+nBTMzKzBScHMzBqcFMzMrMFJwczMGpwUzKZR4ruSTm3a9iZJ15cZl1kRPCTVrAVJLyVZhmIJyUqTm4DXRsT9s9jn3EgWWjOrLCcFsxlIuoJkvaUDgCcj4jJJ5wHvIFmG+Dbgwoh4VtJqkiU4RoAvRcQH0308DPwj8FrgYxFR6npHZp3MLTsAswr7AHAH8DQwntYe3gD8fkQ8kyaCs0mux7EyIn6RLiF9s6RrImJqSY5fRcSyMgpg1isnBbMZRMSvJH0J+PeI+LWkV5FcnGUiWZqGEfYsS3yOpAtI/qcOJ7kgylRS+FKxkZv1z0nBrL1n2XO9DQGfiYi/an6CpKOBi4ATImJ7utpu82Upf1VIpGYZ8Ogjs+59G3iTpEMA0ouxLCRZjPBJYEe6nPEpJcZoNiuuKZh1KSLukvQB4NuS9iG5nOOfkqzi+iOSlUIfAG4tL0qz2fHoIzMza3DzkZmZNTgpmJlZg5OCmZk1OCmYmVmDk4KZmTU4KZiZWYOTgpmZNfx/HheHmFHkKiYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the momentum factor\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.plot(factor, 'o')\n",
    "ax.set_xlabel('Year')\n",
    "ax.set_ylabel('MOM factor')\n",
    "fig.savefig('../output/mom_factor.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regression\n",
    "\n",
    "Here, we run a regression using MKT, SMB, HML factors together with the newly constructed momentum (MOM) factor to estimate the $\\beta$-factors in the following relation:\n",
    "\n",
    "$R_{p,t}^e = \\alpha + \\beta_M R_{M,t}^e + \\beta_{SMB} R_{SMB,t}^e + \\beta_{HML} R_{HML,t}^e + \\beta_{MOM} R_{MOM,t}^e + \\epsilon_t$\n",
    "\n",
    "For the excess stock returns, we'll use Apple excess returns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AAPL_E</th>\n",
       "      <th>MKT</th>\n",
       "      <th>SMB</th>\n",
       "      <th>HML</th>\n",
       "      <th>MOM</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2013-12-31</th>\n",
       "      <td>6.865672</td>\n",
       "      <td>28.05</td>\n",
       "      <td>6.03</td>\n",
       "      <td>1.71</td>\n",
       "      <td>54.163134</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-01-31</th>\n",
       "      <td>23.311377</td>\n",
       "      <td>25.29</td>\n",
       "      <td>5.14</td>\n",
       "      <td>0.45</td>\n",
       "      <td>52.247775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-02-28</th>\n",
       "      <td>14.494308</td>\n",
       "      <td>20.68</td>\n",
       "      <td>6.45</td>\n",
       "      <td>-1.67</td>\n",
       "      <td>51.737800</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-03-31</th>\n",
       "      <td>19.797606</td>\n",
       "      <td>21.30</td>\n",
       "      <td>5.99</td>\n",
       "      <td>-1.78</td>\n",
       "      <td>53.810023</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-04-30</th>\n",
       "      <td>21.746144</td>\n",
       "      <td>20.18</td>\n",
       "      <td>6.52</td>\n",
       "      <td>2.67</td>\n",
       "      <td>52.134021</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               AAPL_E    MKT   SMB   HML        MOM\n",
       "Date                                               \n",
       "2013-12-31   6.865672  28.05  6.03  1.71  54.163134\n",
       "2014-01-31  23.311377  25.29  5.14  0.45  52.247775\n",
       "2014-02-28  14.494308  20.68  6.45 -1.67  51.737800\n",
       "2014-03-31  19.797606  21.30  5.99 -1.78  53.810023\n",
       "2014-04-30  21.746144  20.18  6.52  2.67  52.134021"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get Apple excess returns\n",
    "aapl = (df['AAPL'] - df['RF']).rename('AAPL_E')\n",
    "comb_df = pd.concat([\n",
    "    aapl,\n",
    "    df['MKT'],\n",
    "    df['SMB'],\n",
    "    df['HML'],\n",
    "    factor.rename('MOM')\n",
    "], axis=1)\n",
    "comb_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>AAPL_E</td>      <th>  R-squared:         </th> <td>   0.568</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.542</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   21.73</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Mon, 10 Feb 2020</td> <th>  Prob (F-statistic):</th> <td>1.79e-11</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>13:45:14</td>     <th>  Log-Likelihood:    </th> <td> -286.27</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    71</td>      <th>  AIC:               </th> <td>   582.5</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    66</td>      <th>  BIC:               </th> <td>   593.9</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     4</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   40.9497</td> <td>   17.088</td> <td>    2.396</td> <td> 0.019</td> <td>    6.832</td> <td>   75.067</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>MKT</th>       <td>    1.8038</td> <td>    0.258</td> <td>    6.989</td> <td> 0.000</td> <td>    1.288</td> <td>    2.319</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>SMB</th>       <td>   -0.0784</td> <td>    0.335</td> <td>   -0.234</td> <td> 0.816</td> <td>   -0.747</td> <td>    0.591</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>HML</th>       <td>   -0.8120</td> <td>    0.275</td> <td>   -2.952</td> <td> 0.004</td> <td>   -1.361</td> <td>   -0.263</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>MOM</th>       <td>   -0.8444</td> <td>    0.297</td> <td>   -2.847</td> <td> 0.006</td> <td>   -1.436</td> <td>   -0.252</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 0.080</td> <th>  Durbin-Watson:     </th> <td>   0.630</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.961</td> <th>  Jarque-Bera (JB):  </th> <td>   0.108</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.070</td> <th>  Prob(JB):          </th> <td>   0.947</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.870</td> <th>  Cond. No.          </th> <td>    564.</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                 AAPL_E   R-squared:                       0.568\n",
       "Model:                            OLS   Adj. R-squared:                  0.542\n",
       "Method:                 Least Squares   F-statistic:                     21.73\n",
       "Date:                Mon, 10 Feb 2020   Prob (F-statistic):           1.79e-11\n",
       "Time:                        13:45:14   Log-Likelihood:                -286.27\n",
       "No. Observations:                  71   AIC:                             582.5\n",
       "Df Residuals:                      66   BIC:                             593.9\n",
       "Df Model:                           4                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     40.9497     17.088      2.396      0.019       6.832      75.067\n",
       "MKT            1.8038      0.258      6.989      0.000       1.288       2.319\n",
       "SMB           -0.0784      0.335     -0.234      0.816      -0.747       0.591\n",
       "HML           -0.8120      0.275     -2.952      0.004      -1.361      -0.263\n",
       "MOM           -0.8444      0.297     -2.847      0.006      -1.436      -0.252\n",
       "==============================================================================\n",
       "Omnibus:                        0.080   Durbin-Watson:                   0.630\n",
       "Prob(Omnibus):                  0.961   Jarque-Bera (JB):                0.108\n",
       "Skew:                          -0.070   Prob(JB):                        0.947\n",
       "Kurtosis:                       2.870   Cond. No.                         564.\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Run regression\n",
    "model = sm.ols('AAPL_E ~ MKT + SMB + HML + MOM', data=comb_df).fit()\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZgU1bn48e/LMMCAyMhiFJAlhriBAk6MCi6BqLghaiTuxKvxJt7ELYqQGw16zXVcokl+ifF6MeJNFMWNYNwVNQFjFBB3XEFZRFkcUBhlmHl/f1QN9DRVM9XdVV3V1e/neeaZ6WWqTnVX96lzznveI6qKMcYYE0S7uAtgjDGmdFilYYwxJjCrNIwxxgRmlYYxxpjArNIwxhgTmFUaxhhjArNKwySSiBwqIsviLkdYRGSAiKiItHdvPyoiE4qw3yki8peo92PKh1UaxpOIPCsin4lIx4DPb/GlWIpEZImI1IvIFyLyiYjcLiLbRbEvVT1SVe8IWKbvRlEGt2Juco/3cxF5W0TOyuH/p4nI1VGUzSSXVRpmGyIyADgIUGBsrIUpvmNVdTtgOPAt4BfZTxBHWj47K9zj3R64CPhfEdmtGDsu5QuMcpaWE9+E60zgBWAa0KILRUSqROTXIvKhiKwTkTkiUgX83X1KnXvlekB214hHF81ZIvKWe5X7gYj8e5DCicgtInJD1n1/FZGL3b8vE5HlGVfPo3N9AVR1OfAoMNjd5rMi8isRmQtsBL4uIt1E5DYR+djd39UiUuE+v0JEbhCR1SLyAXB0VnmfFZFzMm7/MOO1eFNEhovIn4F+wEPuazrRfe7+IvK8iNSJyCsicmjGdgaKyHPudp4EegY8XlXVR4C1wN4Z29tdRJ4UkbXuaznevf9c4DRgolu2h9z7VUS+kfH/W1ojzV2O7vuzErg9476ficin7mt5Vsb/H+W+Hp+7r/ElQY7HREhV7cd+WvwA7wHnAfsCDcDXMh77A/As0AeoAA4EOgIDcFom7TOeOwX4S8btFs/B+SLdFRDgEJwv4+HuY4cCy3zKdzCwFBD39g5APdAb2M19rHfGPncNeNxLgO+6f+8CvAH8l3v7WeAjYC+gPVAJzAT+B+gC7Ai8CPy7+/wfAYvc7XQHnsk69meBc9y/TwKW47RsBPgG0D+7TO7tPsAa4Cici77D3Nu93Mf/CdzovicHA59nvgdZx7vlNXa3NRZoAoa593VxX8uz3GMeDqwG9nIfnwZcnbVNBb6RcXvLc9z9bQaudctXlXHfVe5repR7Huzg/s/HwEEZ7/PwuD8f5f5jLQ3TgoiMBPoDM1R1PvA+cKr7WDvg34ALVHW5qjaq6vOq+lU++1LVh1X1fXU8BzyB0y3Wln/gfDk1P/d7wD9VdQXQiPOFtKeIVKrqElV9P4dizRSROmAO8Bzw3xmPTVPVN1R1M05FcCRwoapuUNVPgZuAk93njgd+o6pLVXUtcE0r+zwHuE5VX3Jfi/dU9UOf554OPKKqj6hqk6o+CcwDjhKRfjgVz+Wq+pWq/h14qI3j7e0ebz3wIHCxqr7sPnYMsERVb1fVzaq6ALgf5/XOVxPwS7d89e59DcBVqtqgTmvnC5zKv/mxPUVke1X9zC2DiZFVGibbBOAJVV3t3r6LrV1UPYFOOBVJwUTkSBF5we36qMO5ymyzO0VVFbgbOMW961TgTvex94ALcVo5n4rI3SLSO4dijVPValXtr6rnZXyxgXPV3aw/zpXxx243UR1Oq2NH9/HeWc/3qwTAaY0EfU37Ayc179Pd70hgZ3efn6nqhoD7BWdMoxpnTON3wKisfX07a1+nATsFLKuXVar6ZdZ9a9yKuNlGoDkA4USc8+JDt9vtgAL2bUJglYbZwh2bGA8cIiIr3X7ni4B9RGQfnK6JL3G6lLJ5pUveAHTOuL3ly0acqKz7gRtwur+qgUdwumeCmA58T0T6A992t+UURPUuVW1uMSlOd0gYMo9xKfAV0NOtZKpVdXtV3ct9/GOcyqBZv1a2uxTv1zR7n83P/XPGPqtVtYuq1rr73EFEugTc79adOK3Fy4AhIjIuY1/PZe1rO1X9sU/ZwPnC93zPW/mf1sr1kqoeh1MZzwRm5PL/JnxWaZhM43C6d/YEhro/e+B0B52pqk3An4AbRaS3O9h7gFsBrMLpevh6xvYWAgeLSD8R6QZMznisA0430ipgs4gcCRwetKBuF8oqYCrwuKrWAYjIbiIyyi3TlzjdLo25vhAB9v8xTnfar0VkexFpJyK7isgh7lNmAOeLSF8R2QGY1MrmpgKXiMi+4viGWxkCfELL1/QvwLEicoT7+ndyB5P7ul1a84ArRaSD29V4bA7HtAn4NXCFe9ffgG+KyBkiUun+fEtE9vApGzjv+alu2cbgjFXlxT2G00Skm6o2AOuJ4L00ubFKw2SaANyuqh+p6srmH+D3wGniRD1dArwGvIQTaXMt0E5VNwK/Aua6XRn7u/3t9wCvAvNxvoQAUNXPgfNxvlw/w+limpVjeacD38XpQmvWEajFaRWtxLlC/TmA+wX0Ro77aM2ZOJXfmzjHcB9ONxHA/wKPA68AC4AH/DaiqvfivHZ34Qxcz8QZMwFnLOQX7mt6iaouBY5zj2kVTmvgUrZ+lk/FaXmtBX4J/F+Ox/QnoJ+IHOu+R4fjjNOswHk9mwexAW7DGW+oE5GZ7n0X4FRUzV1ZMynMGcASEVmPE1xweoHbMwVqjj4xxhhj2mQtDWOMMYFZpWGMMSYwqzSMMcYEZpWGMcaYwFKRMKxnz546YMCAuIthjDElZf78+atVtVcu/5OKSmPAgAHMmzcv7mIYY0xJEZG2MgZsw7qnjDHGBGaVhjHGmMCs0jDGGBNYKsY0jDHp0dDQwLJly/jyy+xkuCZfnTp1om/fvlRWVha8Las0jDGJsmzZMrp27cqAAQMQCZr02PhRVdasWcOyZcsYOHBgwduzSqOMzXx5Odc//jYr6urpXV3FpUfsxrhhfeIulilzX375pVUYIRIRevTowapVq0LZnlUaZWrmy8uZ/MBr1Dc4maaX19Uz+YHXAKziMLGzCiNcYb6eNhBepq5//O0tFUaz+oZGrn/87ZhKZIwpBVZplKkVdfU53W9MuXnwwQcRERYtWtTq86ZNm8aKFSvy3s+zzz7LMccck/f/F5tVGmWqd3VVTvcbk1QzX17OiNrZDJz0MCNqZzPz5eWhbHf69OmMHDmSu+++u9XnFVpplBqrNMrUpUfsRlVlRYv7qioruPSI3WIqkTG5ax6bW15Xj7J1bK7QiuOLL75g7ty53HbbbS0qjeuuu44hQ4awzz77MGnSJO677z7mzZvHaaedxtChQ6mvr2fAgAGsXr0agHnz5nHooYcC8OKLL3LggQcybNgwDjzwQN5+uzS7gm0gvEw1D3Zb9JQpZa2NzRVyLs+cOZMxY8bwzW9+k+7du7NgwQI++eQTZs6cyb/+9S86d+7M2rVr6d69O7///e+54YYbqKmpaXWbu+++O3//+99p3749Tz31FD//+c+5//778y5jXKzSKGPjhvWxSsKUtKjG5qZPn86FF14IwMknn8z06dNpamrirLPOonPnzgB07969tU1sY926dUyYMIF3330XEaGhoaGgMsbFKg1jTMnqXV3Fco8KopCxuTVr1jB79mxef/11RITGxkZEhBNPPDFQ6Gr79u1pamoCaDGr/fLLL+c73/kODz74IEuWLNnSbVVqbEzDGFOyohibu++++zjzzDP58MMPWbJkCUuXLmXgwIF0796dP/3pT2zcuBGAtWvXAtC1a1c+//zzLf8/YMAA5s+fD9Ci+2ndunX06eO07KdNm5Z3+eJmlUbEoorsMMY4XazXnDCEPtVVCNCnuoprThhSULfr9OnTOf7441vcd+KJJ7JixQrGjh1LTU0NQ4cO5YYbbgDgBz/4AT/60Y+2DIT/8pe/5IILLuCggw6iomJrhTZx4kQmT57MiBEjaGxsOQ5TSkRV4y5DwWpqajSJizBlz7oG5yqo0JPamDR766232GOPPeIuRup4va4iMl9VWx/Bz2JjGhGKKrLDmNZYTjETJas0ImSzrk2xWU4xEzUb04iQzbo2xWY5xUzUrNKIkM26NsVmrVsTNas0IhRFZIcxrbHWrYmajWlEzGZdm2K69IjdPCP2rHVrwmItjTJic0bSz1q34aioqGDo0KEMHjyYk046acuEvnxkpj6fNWsWtbW1vs+tq6vj5ptv3nJ7xYoVfO9738t731GwlkaKZYZedquqZMOmzTQ0OvNyLKomvcqudfvqDHj6Kli3DLr1hdFXwN7jC9pkVVUVCxcuBOC0007jlltu4eKLL97yuKqiqrRrl9t199ixYxk7dqzv482VxnnnnQdA7969ue+++/I4guhYSyOlslNG19U3bKkwmllUjSl5r86Ah86HdUsBdX4/dL5zf0gOOugg3nvvPZYsWcIee+zBeeedx/Dhw1m6dClPPPEEBxxwAMOHD+ekk07iiy++AOCxxx5j9913Z+TIkTzwwANbtjVt2jR+8pOfAPDJJ59w/PHHs88++7DPPvvw/PPPM2nSJN5//32GDh3KpZdeypIlSxg8eDDg5LE666yzGDJkCMOGDeOZZ57Zss0TTjiBMWPGMGjQICZOnBjasXuxSiOlvEIvvVhUTWmxLsYsT18FDVnncEO9c38INm/ezKOPPsqQIUMAePvttznzzDN5+eWX6dKlC1dffTVPPfUUCxYsoKamhhtvvJEvv/ySH/7whzz00EP84x//YOXKlZ7bPv/88znkkEN45ZVXWLBgAXvttRe1tbXsuuuuLFy4kOuvv77F8//whz8A8NprrzF9+nQmTJiwJSHiwoULueeee3jttde45557WLp0aSjH76V8K41XZ8BNg2FKtfM7xCuTJAhaGVhUTemIasGhkrZuWW73B1RfX8/QoUOpqamhX79+nH322QD079+f/fffH4AXXniBN998kxEjRjB06FDuuOMOPvzwQxYtWsTAgQMZNGgQIsLpp5/uuY/Zs2fz4x//GHDGULp169ZqmebMmcMZZ5wBOGtz9O/fn3feeQeA0aNH061bNzp16sSee+7Jhx9+WNDxt6Y8xzSam7TNVyjNTVoouC80KfxSRmeyqJrSYmlpPHTr63ZNedxfgMwxjUxdunTZ8reqcthhhzF9+vQWz1m4cGGgFOq5ai1PYMeOHbf8XVFRwebNm0Pff7PybGlE3KQtVBhdEF4TCyvbCTt0rrSomhIV98S9RHaNjb4CKrNay5VVzv0R23///Zk7dy7vvfceABs3buSdd95h9913Z/Hixbz//vsA21QqzUaPHs0f//hHABobG1m/fv02adYzHXzwwdx5550AvPPOO3z00UfstlvxL/rKs9KIqEkbhrC6ILxCL68/aR9evuJwFtcezdxJo6zCKDFxTtxLbNfY3uPh2N9Bt10AcX4f+7ui9Bj06tWLadOmccopp7D33nuz//77s2jRIjp16sStt97K0UcfzciRI+nfv7/n///2t7/lmWeeYciQIey777688cYb9OjRgxEjRjB48GAuvfTSFs8/77zzaGxsZMiQIXz/+99n2rRpLVoYxVKeqdFvGuzTpN0FLno9vILlYUTtbM9upT7VVcydNCqGEpmkiDPVfjHPS0uNHo2wUqOXZ0sjxiZtW+LugjDJFefEPTsvTbPyHAhvbrpmTwgCtxUS3iShXBWy5rGto5B+cU3ci2ItblOayrOlAU5lcNHrMKVua5dUxJOEgsg3M25i+5xNKhQ7Y3Maus2TJMzXs3wrjWwJiajKtwsirHUUEhkhY2JXzK6xTp06sWbNGqs4QqKqrFmzhk6dOoWyvdi7p0SkApgHLFfVY0RkIHA30B1YAJyhqpsiL0iCIqry6YIIo8/ZVn0zrck+L5svMMLuDu3bty/Lli1j1apVBW/LODp16kTfvoXNXWkWe6UBXAC8BWzv3r4WuElV7xaRW4CzgT9GXoqIJgkVSxh9zjZ5zAQV5QVGZWUlAwcOLLiMJhqxdk+JSF/gaGCqe1uAUUBzWsc7gHFFKUySIqrySHESRp+zRciYoGxZ2fIVd0vjN8BEoKt7uwdQp6rNc+CXAcW5xPWLqCo0eirXtM15pjhpvrorJHrKImRMUEm/wLBIwujEVmmIyDHAp6o6X0QObb7b46meo2Eici5wLkC/fv3CKdTe48MNsfWpAF5a8hkXvjnI+4RubUC+jbIVGo5pq76Vr1y/ZJN8gWFjc9GKs3tqBDBWRJbgDHyPwml5VItIc2XWF1jh9c+qequq1qhqTa9evYpR3tz5VAC951/nHxob44C8rfpWnvIJ1y52CG4urOssWrG1NFR1MjAZwG1pXKKqp4nIvcD3cCqSCcBf4ypjwXy+6HdmTYvbLQabYx6QL7tV31Kg0K6YfAIgwugOjUrSu85KXdxjGl4uA+4WkauBl4HbYi5P/nwqgBXao8Xtse3mMHHjDJiyBqp2gIoO0JgRZZyQFCcmecLoisn3SzapFxhJ7jpLg0RM7lPVZ1X1GPfvD1R1P1X9hqqepKpfxV2+vHlEZNXTkaebhjKnw/l80PFU5nc4lxsqb6Vvu9WAQv1aUIWq7hQ7a6cpPWF0xcSZPTcKSe46S4MktjTSwyMia0X1CMYveYAqcVoSPeSLbf+vqQE6dIHLFhexsKYUhdEVk7YAiCR3naWBVRpRy4rI2vWmwSABJrgnYG0Pk3xhdMWE+iWba4h5RJLadZYGVmm0JooPQNDKIKEz0S3+PVnCaiWE8iXrFWI+8zx49DKo/yzWSsSExyoNP1F9APyiozIldODb4t+TJ1FdMV4h5k0NzjgdBJ6oapKtPFfuC8Jvdb9MlVXbDlK31TrJrowA2lVCx66JvxqzVQXTI5IW45RqfObitpSAFTKNI5+V+6yl4SdIN1L2TO0gKUCiSldSBBb/ng6RtRiDtKLBxutKnFUafvL5AARNARJ2upIIZV6RthOh0aNlWqqhmeUiu1WxcdPmSLIZv7TrTxk8/xdbIgN9JXS8zgRjlYaf0Vds243kJfMDkGcKkKQOLmdfkXpVGKUcmlkOvFoVfnJuMWZ1xT684UT+3HAOE9vPoLes4TPtQlf5kg6yeev/JHS8zgSXiMl9ibT3eGe8otsugDiT7So6tHxO9gfA7wqqlSurJC/T6jVxDKBCxHJTlQi/99BLTi3G5q7YjOWRJzbcDMDITb/j61/dyb6bbuWShnNZ1tQTm6iaHtbSaE12N1Jbg9xerZM2rqzyXfioGK0TvyvPJlUW1x4d6r5MNIK2HnJuMXp0xXaWTUxsP4NZm0ZuuW9W00jmdz7MAiVSxCqNXLQ1FpHHIHfQweXMSqJbVSUbNm2modHpLvIdyCxwnonl8Cl9fu9hdVUlXTq2z/+iw6fLtbe0TMZp3ZfpY5VG2HIc5A7yxZzdL11X37DN87dpnQSJ5GqjUklbeoly5PceThm7V2EtWZ9AkS8770Sfqqr4x+cSMjM9jazSiFmQL+ag/dItWidtRXIFqFQSNXHM5CWf9zBQSK5PV2znI69i7t4xd0XlufqlCcYqjUIVeEUT5EO9oq7eSZ/efga9ZTUrtCfXbR7PrKaRLbbVotuorUiugOHBlsOn9OX6HgYaZ0vyfKMCVr80bbNKoxAhXdG09aGesN2LTGyYSmc3/r2vrKa2cio0sKXi2KbbqK3FnGJcIdAkW+BJnEmdb2TndqQs5LYQrV3RhGhi5T1bKoxmnWUTkzrM8A999VjLo0UkVx7hwaY8lPz6GnZuR8oqjUIU6Yqmc/1Kz/t7s4bFtUczd9KobVsq2fNMsmPk26pUTNkq+UWM7NyOlHVPFcKnC2glPTlg0sPhDRz7dTVJOydJnF9/cmvdByH2SSd1RrvJz7hhfeiz9G/ssuB6dtRVfCq9WDr8Ur41bEzcRQvG59ye2TiC62tn23laIMtyWwiPjLX12oHLGs5pMdZQ8Kxpr8y42bwy7oahjYH+7EgbCOmYTXy8zreozq8isfPUWz5Zbq17qhBZXUAr6dWiwoDc12sOsh+kYtvnRDCW4pUqgofOd+53hbFGtUmYIo3VFZOdp+Gx7qlCZXQBHTDpYc/VBEJJHZ7Z1TSl2vs5YUeHBAhdtHTpKZTC6CO/sPWH6ka2/c+mBWtphKhoUSfFig4J8OVR8pE2ZlspjD6asN2L1FZOpW+71bQT6NvOCVufsN2LcRet5FilEaKiRZ34RYcMOtxZcXBKtfM7oxspLwG+PEo+0sZsK8roo1dnhHuOBuQXtj6x8p6i7D9NrHsqREVLu+EVHTLocHjlrnBTJwTI2mupRlLIL7IO3GWQ84y2y3MybBjReX5h6373G38WPZUWfmua57oec3a01KDD4d0nkpcqwhRXGBFVeZyjoUU9hfX5SBmLnipnYQxeekVLvXKXU1FMqXM+XFZhlKcwIqryOEcDRz211e1lE/5C02alISIXiMj24rhNRBaIyOHFKJzJQRiDlykMtTQhCeOiJI9zNFB0XoDQ8DYzJJjAgrQ0/k1V1wOHA72As4DaSEtlchfGlVQKQy0jU6QB3ZkvL2dE7WwGTnqYEbWzvZcBLkZZwrgoyeMcDRSdF/RiZ+/xTmvZWs0FCTIQLu7vo4DbVfUVEZHW/sHEIIy0IG1lxjWOiNZryB7w/c7uvbh//vLW17Uo1toReSxlvI08ztFAC4HZxU5RtTkQLiK3A32AgcA+QAXwrKruG33xgrGB8JCkMH1EJCIYVPUa8BXYZrLo2HZz+HmHe9mJ1c6X7qYNUL821LL4imk1vDajp/J8PyxnWn4D4UFaGmcDQ4EPVHWjiPTA6aJKnbI/iZK8sE6SRHBl6zXg61Vh1FZOpTPufAOvL8oQyuIrpvUz2lxEKo9WUKDVCY2nNisNVW0SkU+APUUktfM67CRyJXVhnSSJoBsvSNqVie1nbDNBzVc5dSnmcbETaHVC46nNSkBErgW+D7wJNL/KCvw9wnIVnZ1EJrAw+vez9K6uYrlHxZHZRdVbVgfbWDmGkuZ4sWM50/IXpOUwDthNVb+KujBxspPIBBZBN57fgO+J+/bhmUWrWFFXz6fSi51Yte0/V3WHDl1KoksxjC7gMLbhV0lbzrS2Bak0PgAqgVRXGvmeRGU/DlKuQu7GC5SO5dUN3i2cI69NbCWRKYwu4LC6kQNFZRlPQSqNjcBCEXmajIpDVc+PrFQxyOckKstxkJgiaMpBmwO+JR6oEEYX8PWPv81hjc8xsUPLFOfXP94hp89cUXOmpewzE6TSmOX+pFo+J1HZjYP4zQn46AXLT1UsYbVwYvgiC6MLuGb9k1xTOXVLQEBfcVKcT14PMCqn8rRZSYehWPNoiqjVSkNEKoDDVPX0IpUnVrmeRGU3DuI383ben9gyXJuCD0XqFfGLLLP7tp0IjR7zwnIZR5jc4d6tIceuzrKJyR3uBa4ptLjhC7CQWalpNY2IqjYCvUSkQ9g7FpFdROQZEXlLRN4QkQvc+7uLyJMi8q77e4ew9x2WsluAyDf2P+uLwPJVJVuRcow1d98ur6tHwbPCyHUc4Wt4R5D53R+7FM5WD9I9tQSYKyKzgA3Nd6rqjQXuezPwM1VdICJdgfki8iTwA+BpVa0VkUnAJOCyAvcVirbSPEDKB9P85id4KeEPRbHEFkRRpC8yr+5bgAoRmlTpXV3Fb/Z8l289ewn8NVg3mficg5KUeSnZ3X5VO/jM2E9IefMQJGHhCuBv7nO7ZvwURFU/VtUF7t+fA2/hpCs5DrjDfdodOCG/scu+alpeV8/985dz4r596FNdhQB9qqtyz/NfSrwSzuGThqyEPxTF4HU+zXnwZjZeu3v0q9oVaTlXv27aJlUW1x7N3KNW863Xftl6dtpsSU5x7pVt96vPoSKroyYp5c1TkBnhV0ZdCBEZAAwD/gV8TVU/dvf9sYjs6PM/5wLnAvTr1y/qIvoOej+zaBVzJ+U2AFeygqwYCCX/oSiG7PNpbLs5XCVT6VyfkSIkqrGhCCYnemkzjD2f/v4kR5B5HU9TQ0nNowkiyIzwZ9g2DQ6qGso3pYhsB9wPXKiq64Mm0FXVW4FbwUlYGEZZWlN2g95+vKJ3+u2fzA9xgmWfN54pQhrq4dHLwn9ti/TF22YYe77dZElNdeNX7vrP4LLFxS1LhIKMaVyS8Xcn4ESc8YiCiUglToVxp6o+4N79iYjs7LYydgY+DWNfhbIZpK3I+hDPfHk519fOtgmPrcg+n3xThNSv3donHrT1ESScNvuLt3lNjhArkTbD2IPm8CqVeQ5lsrRAm2Maqjo/42euql4MfLvQHbtrctwGvJU1qD4LmOD+PQH4a6H7CsOlR+xGVWVFi/tSPeidJ6+++skPvOa9eFAZyz6fVmjPYP/YVpRTkFXswvifgMYN68PcSaOcMYxJo1pePAQZn4iwbKFL8nhLiIIs99o946eniBwB7BTCvkcAZwCjRGSh+3MUzqqAh4nIu8BhJGSVwHHD+nDNCUPKZ9A7T4HXdC5z2efT1A6ns7miU7B/bq37Jp9w2riW+Q2yBGspLUFcJkvKBumemo8zpiE43VKLcdbYKIiqzsE39IbRhW4/CkWZQVribOwnuJbn09Hw6l4tu2F8F1hqpbsjn3GCOOcStDU+UWrzHJI63hKiIJXGHqr6ZeYdItIxovKYEmdjPwXwGmfINcopn371GPvi25yrkvRxglIZbwlRkHkaz3vc98+wC2LSIbVjP80DxVHPociUT3dHPv3qMfXFBxr/SvI4QSmNt4TIt6UhIjvhTLarEpFhbO1K2h7oXISymRJU1OyhxRJn0rlcuzvyCaeNae5DoISfpTYvo8TzSgUh6pEPBkBEJuCk9KgB5mU8tB64IyNENnY1NTU6b968tp9oTD5uGuzTRbILXPR6/tstw66NTAMnPbztBDDguHZz+G2vh3J6XWJJyTKlGo8pbIDAlLpo9x0SEZmvqjW5/I9vS0NV7wDuEJETVfX+gktXgmyBJQNEMxibwpTZufIa/xrbbg61HW6Dde7SPQFel9jWtUn6eEtEgoxpzBWR20TkUQAR2VNECo6eSrpUzjeIo18+DaLI1VRKoaQR8aZNvFMAABa+SURBVBr/uqxyBlXZi4S28brEFuad5PGWCAWpNG4HHgd6u7ffAS6MrEQJkbr5BmU6aBeKoJPQcqmQSy2UNAJec596yxrvJ7fyusQW5u0RqPDSkCsZ8UhPBk56mBG1s0v7ItNHkJDbnqo6Q0QmA6jqZhHZNt9xyqRuvkExB+08+upnNo4o3a6+tgZj8+lqKmbXRoLHTraZ+3RT7q9LrGHeGYEKW7vJnLKkdfnnIC2NDSLSA3fER0T2B9ZFWqoESN0CS8W6svVo0Wz+60+Z8+DNpd3Vt/d4Z9B7Sp3zu9BZy8Xq2ii1FmYer0tSwrxT1zvhI0ilcTFOPqhdRWQu8H/A+ZGWKgGSciKGpkhrKHh9gbZv/JILubvFfan6MOVTIRcr5USpjZ3k8bokJcVP6nonfARZT2OBiBwC7IYzV+NtVW2IvGQxS918gyKtoeD3RenVV52aD1O+XU3FSDlRimMnebwuSUjxUy7ZEIKMaaCqm4E3AETkMBGZqKqHRVqyBEjCiRiaYk2S8vkCXaE9trkvNR+mYlXI+SjTsNA4tLl+SNhiGqtqbUb4KOAWnKipmcB/43RNCfCryEtmwleMK1uPL9DNFZ34TdPJLZ5W0l192ZI8aznJFVrKFLV3IsZ5Pq3NCH8ZuAgnz9SROBXG5ar620hLlAebEZ4waYueKnUJjp4yeQopS0E+M8JbqzQWqOrwjNvvq+quuWy8WKzSMMaUtFwr9pBSmISaRgSoFpETWm5/6+0k5Z4yxpiSlfR5PllaC7l9Djg24yfz9jGRl8wYY8pBkuf5eGgtYeFZke/dGGPKXb7zfCBZ0VPGxM2yDJuykOR5Ph6CzAg3puhSmWXYBDLz5eWMqJ1d/KR/cWWBLrFsudbSMIkUaFU3kzpFXRsjM2KpagfY9AU0bnIeK/bqjFAyYdFtVhoichLwmKp+LiK/AIYDV6vqgshLZ8pWueTxMS0V7WIhO2Kpfu22zynm0q0xdTXlI0hL43JVvVdERgJHADcAfwS+HWnJImb95cmWijw+WbH3L+36Uy58c5Cdc60o2sWCV8SSlyTn6IpJkDGN5mr/aOCPqvpXoEN0RYqe9ZcnX8lnGfZIST54/i/Yd/2Tds61omhLEgStDEKY9xDbGE1EglQay0Xkf4DxwCMi0jHg/yVWueS9L2VJSXedN48r2SrZxI2Vt/BBx1OZ0+F8Dmt8zs65LEW7WAhSGYQwGJ3GC9QgX/7jcZZ7HaOqdUB34NJISxUx6y8vDeMq5jK34/ks7nQaczuez7jlvy6dNc59rmTbSxPtBPq2W01t5VRq1j9Z5IIlW9EuFjwilhqlPXV0pUmFlfTipSFXFjzOkMYL1CBjGjsDD6vqVyJyKLA3TvLCkpWK/vK080qtMO+2rY8XM7olH36x9xk6yyYmd7gXuKY4ZSoRRVmSICtiaWPVTlyx4UTu23TglqdUvVTBNbssL6gsabxADdLSuB9oFJFvALcBA4G7Ii1VxEq+v7wU5RoDH2SgMskr0HnF3nv4GquLUBjjKWMJ38P05hYVBoTTIkjdstEEqzSa3EWYTgB+o6oX4bQ+SlbJ95eXmnzWqQ46UJnU6JbsZUulwvNpYoshJUJULYI0XqAG6Z5qEJFTgDNxkhUCVEZXpOJI1ap8SeeTkG3lAz/ngLu6eIefBuje2fK8pMqMvc/uboNEz/otN1F1Wee9MFOC10AJUmmcBfwI+JWqLhaRgcBfoi2WSRWf1sCOurpFRAlkzPr1WnEuWyl96ZbYrN9yE+VSrUEuUDPnjU3Y7kV+obfQvvFL58GEjd/5LsLU4kkiVUA/VU3kkL8twpRwPquMLWvqychNvwNgbLs5/LzDvezE6q1fqNDyS3bQ4fDuE/alayIRaMJvBBM2s1OnzOlwPn3beYx15bgqXxBhL8LUvNFjcWaBdwAGishQ4CpVHZtfMU3Z8Wg1bNQOXLfZ+cIf224OtZVT6UxW3p9jfxf6h8QYP222CDwi+gbP/wX7NpzDckbmnScrOyy3t/gERyRk/C7IQPgUYD+gDkBVF+JEUBkTTNag8Ep6ManhHGY1jQRgYvsZdJZNLf8nyZFRpjz5TNic2H5rQEc+EVfZg+0rtKf3E7PG7+KaaR6k0tisquuy7mu7T8uYTBnhjS8c9xxPVhyy5aGkX1kZA/iej71lTYvbuUZcZQ+2X7d5PBs1K1NT1vhdnDPNg1Qar4vIqUCFiAwSkf8HPB9xuUyKZYc8rxTvK6uNVTsVt2DGtMYnUm+F9mhxO9eIq+yw3FlNI7lCz2Vj1c6AOC30Y3/XYvwuzpnmQSqNnwJ7AV/hTOpbB1wYZaFM+o0b1oe5k0axuPZobm1/+jZXVhu1A9c1fD+m0hnjwWPCZn3G2BzkF3HlNW9s5PHn0fmyRTClzmmhZwV8xDnTvM2BcFXdCPyn+1M0IjIG+C1QAUxV1dpi7t8Uzx1f7Mfadk7fcG9ZwwrtwXWbx/PQV/sxJe7CGdPMI2z69V1/yvw3ByEFprvPdd5YnKmQgkRPPQmc5CYrRER2AO5W1SOiKpSIVAB/AA4DlgEvicgsVX0zqn2a+PSurmJW3UhmbRrZ4v4+JZxqwaRU1mJJ3wLmxhBHGuW8krYE6Z7q2VxhAKjqZ8CO0RUJcKK13lPVD1R1E3A3cFzE+zQxSWOqBWOiFGcqpCAzwptEpJ+qfgQgIv2JPnqqD5A5G2wZJb5SoPGXd6oFY8pYXKmQglQa/wnMEZHn3NsHA+dGVyQAxOO+FhWViJzbXI5+/fpFXBwTNcsFZkxpCDIQ/piIDAf2x/kyv0hVo87nvAzYJeN2X2BFVrluBW4FJ41IxOUxpgVbY96UqzbHNETkbFVdrap/U9WHgM9E5JcRl+slYJCIDBSRDsDJwKyI92lMIGlcwtOYoIIMhI8WkUdEZGcRGQy8AHSNslDu+h0/wVlm9i1ghqq+EeU+TRnLcYGooBOr4krzYBIu1wXJEiZI99SpIvJ94DVgI3CKqs6NumCq+gjwSNT7MWXOa1nZNtJQB5lYlZ25NN9kdiZl8jjfkiZI99Qg4AKcZV+XAGeISOeIy2VMcfgsENVaskS/CVTtRLa0Kq586I3Y0jyYBMvjfEuaIN1TDwGXq+q/A4cA7+KMORhT+vySIraSLNFrXglAo+qWMY7PNjZ4/m8x0jyYBMvjfEuaIJXGfqr6NIA6fg2Mi7ZYxhSJ33KxrSwjmz2xqkK8IsS9FSPNg0mwPM63pPGtNERkIoCqrheRk7IePivSUhlTLB5J6IIsI5uZcLEpwOqXYLPcDXmfb0nSWkvj5Iy/J2c9NiaCshhTfFkLRHmloW6LX+uhuqoyljQPJsFCON/i1lr0lPj87XXbmNKVlYQuV37J46aM3csqCbOtAs+3uLVWaajP3163jUmmV2e0SGXN6Cvy+sC2NgPccmelVxQz/0s9m4CoT3+siDQCG3BaFVU4czRwb3dS1cqilDCAmpoanTdvXtzFMEmTHRMP0K4SOnaF+s8CVyLZcy7AaUlYd1O6RfG+J+1cEpH5qlqTy//4jmmoaoWqbq+qXVW1vft38+3EVBjG+PKKiW9qgPq1gG6dWBXSDHCTLlG872k4l4KE3BoTvmKkUggS+x5gYlWcS2ua+Pi9v8vr6vNODRP0XEpyChqrNEzxNXcbrVtKLlf8OQsa+95G5eIXHWVzLtKttfc330SVQc6lpCfEtErDFF+xUil4xcR7aaNysZUFy5PfzP9MuXYtBTmX/LqwfjbjlUS0PIIswmRMuIqVSqF5gLs5eqpqB9j0BTRu2vqcgBP5wKKjyk32++4XMppLN2WQc8lve41u0FLcyS99o6dKiUVPlZibBrtdU1m67QIXve77b6GEKoYUgmvKz4ja2Sz3+ELvU13F3EmjIt9PFPsNNXrKmMjkkUohtH7evcc7FdOUOue3VRgmoGJ1UwbpFoP4AjGs0jDFl0cqhTSEKprSlp2oMqrUMEETYsYViGFjGiYeOaZSKGrYq3VhGR/jhvUpyjhC5n78JgTGFYhhlYYpCb2rqzz7eUO/2krBymomXZIWiGGVhikJfkkBQ7/aai0c2CoNE5NitXCCsErDlISiXW2lYGU1Y6JklYYpGflcbeUcptutr084cOmsrGZMlCx6yqRWXmG6KVhZzaRQMXK1BWSVhkmtvMJ0U7CymkmZYuVqC8i6p0xq5R2mW+Irq5mUSVhwhrU0TGpZdlqTCgkLzrBKw6SWZac1qeAXhBFTcIZVGia1ipX2wZhIJSw4w8Y0TKolaVKUMXnJTvEfc2obqzSMMSbpEhScYd1TxhhjArNKwxhjTGBWaRhjjAnMKg1jjDGBWaVh0iVBOXqMSSOLnjLpYQsoGRM5a2mY9GgtR48xJhRWaZj0SFiOHmPSyCoNkx4Jy9FjTBrFUmmIyPUiskhEXhWRB0WkOuOxySLynoi8LSJHxFE+U6ISlqPHmDSKq6XxJDBYVfcG3gEmA4jInsDJwF7AGOBmEanw3YoxmWwBJWMiF0v0lKo+kXHzBeB77t/HAXer6lfAYhF5D9gP+GeRi2hKVR45enJeR9yYMpaEMY1/Ax51/+4DLM14bJl73zZE5FwRmSci81atWhVxEU1a5bWOuDFlLLJKQ0SeEpHXPX6Oy3jOfwKbgTub7/LYlHptX1VvVdUaVa3p1atX+Adgki+EiXx+64j/bMYrDJz0MCNqZ1sFYkyGyLqnVPW7rT0uIhOAY4DRqtpcMSwDdsl4Wl9gRTQlNCUtpIl8fuuFN7qnZHPLA7AuK2OIL3pqDHAZMFZVN2Y8NAs4WUQ6ishAYBDwYhxlNAkX0kS+IOuF1zc0cv3jb+e0XWPSKq4xjd8DXYEnRWShiNwCoKpvADOAN4HHgP9Q1Ub/zZiyFdJEPq91xL34tUiMKTdxRU99o5XHfgX8qojFMaWoW1+nSyrLxqqdOKx2duBIqObHmqOn2ols6ZrKFKRFYkw5SEL0lDG585jIt7miE1dsODHnSKhxw/owd9IoFtceza/H77NNy6OqsoJLj9gt7CMwpiRZpWFKk8dEvqvlR9y36cAWT8t1PGLcsD5cc8IQ+lRXIUCf6iquOWGIDYKXK0u1vw1Rj6Z4qampqdF58+bFXQwTs4GTHvaMzxZgce3RxS6OKXXZEXrgtG5TlGVAROarak0u/2MtDZMafuMONh5h8mKp9j1ZpWFSwysSysYjTN4s1b4nW7nPpEZ2JJTlkTIF8YnQK/dU+1ZpmJLll2jQKgkTitFXeI9plHmqfas0TElqTjTYnDfK0n2Y0DUPdj99ldMl1a2vU2GkZBA8X1ZpmJLkl2jw+sfftkrDhCePVPtpZwPhpiT5pfWwdB/GRMsqDVOSLLzWmHhYpWFKkoXXGhMPG9MwJcnCa42Jh1UapmRZeK0xxWfdU8YYYwKzSsMYY0xgVmkYY4wJzCoNY4wxgdlAuCkrfvmqjDHBWKVhyoblqzKmcNY9ZcpGa/mqjDHBWKVhyoblqzKmcFZpmLJh+aqMKZxVGqZsWL4qYwpnA+GmbFi+KmMKZ5WGKSuWr8qYwlj3lDHGmMCs0jDGGBOYVRrGGGMCs0rDGGNMYFZpGGOMCUxUNe4yFExEVgEfFrCJnsDqkIqTBHY8yWbHk2xpOx7wP6b+qtorlw2lotIolIjMU9WauMsRFjueZLPjSba0HQ+Ee0zWPWWMMSYwqzSMMcYEZpWG49a4CxAyO55ks+NJtrQdD4R4TDamYYwxJjBraRhjjAnMKg1jjDGBpbLSEJFdROQZEXlLRN4QkQvc+7uLyJMi8q77ewf3/t1F5J8i8pWIXOKxvQoReVlE/lbsY3H3H9rxiMgSEXlNRBaKyLwUHE+1iNwnIovc7R1QqscjIru570vzz3oRubBUj8d97CJ3G6+LyHQR6VTs44ngmC5wj+eNON6fPI/nNBF51f15XkT2ydjWGBF5W0TeE5FJbe5cVVP3A+wMDHf/7gq8A+wJXAdMcu+fBFzr/r0j8C3gV8AlHtu7GLgL+FupHw+wBOiZlvcHuAM4x/27A1BdyseTsc0KYCXO5KuSPB6gD7AYqHJvzwB+UMrnHDAYeB3ojLO0xFPAoBI4ngOBHdy/jwT+lXGevQ983f38vALs2dq+U9nSUNWPVXWB+/fnwFs4J/BxOF8yuL/Huc/5VFVfAhqytyUifYGjgalFKLqnMI8nCcI6HhHZHjgYuM193iZVrSvKQWSI6P0ZDbyvqoVkOshLyMfTHqgSkfY4X7QrIi6+pxCPaQ/gBVXdqKqbgeeA44twCC3kcTzPq+pn7v0vAH3dv/cD3lPVD1R1E3C3uw1fqaw0MonIAGAY8C/ga6r6MTgvOs7VRFt+A0wEmiIqYk5COB4FnhCR+SJyblTlDKrA4/k6sAq43e0+nCoiXSIsbptCeH+anQxMD7t8uSrkeFR1OXAD8BHwMbBOVZ+IsrxBFPgevQ4cLCI9RKQzcBSwS3SlbVsex3M28Kj7dx9gacZjy9z7fKW60hCR7YD7gQtVdX0e/38M8Kmqzg+9cHko9HhcI1R1OE4T9T9E5ODQCpijEI6nPTAc+KOqDgM24DTJYxHS+4OIdADGAveGVbY8y1Ho52cHnKvWgUBvoIuInB5uKXMuU0HHpKpvAdcCTwKP4XTnbA61kDnI9XhE5Ds4lcZlzXd5PK3VeRiprTREpBLnxbxTVR9w7/5ERHZ2H98Z+LSNzYwAxorIEpxm2ygR+UtERW5VSMeDqq5wf38KPIjTPC26kI5nGbBMVf/l3r4PpxIpurDeH9eRwAJV/ST8kgYT0vF8F1isqqtUtQF4AKdvPRYhfoZuU9XhqnowsBZ4N6oytybX4xGRvXG62Y9T1TXu3cto2VLqSxtdiKmsNEREcPq531LVGzMemgVMcP+eAPy1te2o6mRV7auqA3C6C2aratGvlMI6HhHpIiJdm/8GDsdpbhdViO/PSmCpiOzm3jUaeDPk4rYprOPJcAoxdk2FeDwfAfuLSGd3m6Nx+t6LLsz3SER2dH/3A04ghvcq1+Nxy/oAcIaqvpPx/JeAQSIy0G3hnuxuw1+ho/hJ/AFG4jSxXgUWuj9HAT2Ap3GuDJ4GurvP3wmnxl0P1Ll/b5+1zUOJL3oqlOPBGQN4xf15A/jPUj4e97GhwDx3WzNxI0RK+Hg6A2uAbmn4/ABXAotwLk7+DHRMwTH9A+fi5BVgdIkcz1Tgs4znzsvY1lE40VfvB/lOsDQixhhjAktl95QxxphoWKVhjDEmMKs0jDHGBGaVhjHGmMCs0jDGGBOYVRrG5Egcc0TkyIz7xovIY3GWy5hisJBbY/IgIoNx0nwMw8kUuhAYo6rvF7DN9uokwTMmsazSMCZPInIdTr6rLsDnqvpfIjIB+A+cNNPPAz9R1SYRuRUnxUkVcI+qXuVuYxnwP8AY4DeqGmu+KWPa0j7uAhhTwq4EFgCbgBq39XE8cKCqbnYripNx1mKZpKpr3RThz4jIfaranPJkg6qOiOMAjMmVVRrG5ElVN4jIPcAXqvqViHwXZ+GeeU5qIKrYmnb6FBE5G+cz1xtnwZzmSuOe4pbcmPxZpWFMYZrYutaKAH9S1csznyAig4ALgP1Utc7NlJy57OmGopTUmBBY9JQx4XkKGC8iPQHchXr64SSL/BxY76arPiLGMhpTEGtpGBMSVX1NRK4EnhKRdjhLhf4IJwvvmziZXj8A5sZXSmMKY9FTxhhjArPuKWOMMYFZpWGMMSYwqzSMMcYEZpWGMcaYwKzSMMYYE5hVGsYYYwKzSsMYY0xg/x/OvNhrrnwYDwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get the parameters, plot the results\n",
    "param_set = model.params.values\n",
    "\n",
    "def prediction(param_set, df):\n",
    "    return param_set[0] + param_set[1]*df['MKT'] + param_set[2]*df['SMB'] + param_set[3]*df['HML'] + param_set[4]*df['MOM']\n",
    "\n",
    "pred = prediction(param_set, df=comb_df)\n",
    "actual = comb_df['AAPL_E']\n",
    "\n",
    "err = (pred-actual)*100/actual\n",
    "\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.plot(actual,'o',label='Actual')\n",
    "ax.plot(pred,'o', label='Prediction')\n",
    "ax.set_xlabel('Year')\n",
    "ax.set_ylabel('Excess Returns')\n",
    "ax.set_title('Actual vs. Predicted Returns')\n",
    "ax.legend()\n",
    "fig.savefig('../output/pred_hw3q3.pdf')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([40.94965675,  1.80380612, -0.07836873, -0.81202584, -0.84441395])"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
